This sets up an agent-based model that exhibits complexity. We will use it to for the following workshop, where you will attempt to optimize over the output from the model.
The model is a square grid of cells, wrapping on both axes, with one agent per cell. Each agent has a value that it can adjust up or down. An agent wants to have a higher value than three of its neighbors, but does not want to be the highest. Agents will:
# From https://stackoverflow.com/questions/2453326/fastest-way-to-find-second-third-highest-lowest-value-in-vector-or-column
maxN <- function(x, N=2){
len <- length(x)
if(N>len){
warning('N greater than length(x). Setting N=length(x)')
N <- length(x)
}
sort(x,partial=len-N+1)[len-N+1]
}
# End sourced code
# Gets the change the agent performs based on its neighbors (not the new value)
getAgentChange <- function(agentValue, neighborValues)
{
neighborValues = c(neighborValues, agentValue)
if (agentValue >= max(neighborValues))
return (-3)
if (agentValue < mean(neighborValues))
return (1)
if (agentValue >= mean(neighborValues) && agentValue < maxN(neighborValues, 2))
return (1)
return (0)
}
# Builds a world with the specified height and width.
# startingValues should be a vector, with one value
# per cell (applied down each column, then across each row).
buildWorld <- function(height, width, startingValues = NULL)
{
if (is.null(startingValues))
{
startingValues = rep(100, height * width)
}
densityValue = 1 / (height * width)
sumStartValues = sum(startingValues)
world <- data.frame(X=0, Y=0, Value=startingValues[1], Density = densityValue)
i = 1
for (x in 1:width)
{
for (y in 1:height)
{
world[i,]$X = x
world[i,]$Y = y
world[i,]$Value = startingValues[i]
world[i,]$Density = startingValues[i] / sumStartValues
i = i + 1
}
}
return (world)
}
# Create and graph an initial world
world = buildWorld(16, 16, seq(100, 125.6, 0.1))
ggplot(world, aes(X, Y, z = Density)) +
geom_raster(aes(fill = Density)) +
geom_contour(colour = "white")
# Updates the passed in world according to the updateFunction
# Neighbors are only the adjacent four. If includeDiagonals is
# true, then all eight adjacent cells are counted as neighbors.
updateWorld <- function(world, updateFunction, includeDiagonals = FALSE)
{
newWorld = buildWorld(max(world$X), max(world$Y))
for (x in 1:max(world$X))
{
xLower = x - 1
if (xLower == 0) xLower = max(world$X)
xUpper = x + 1
if (xUpper > max(world$X)) xUpper = 1
for (y in 1:max(world$Y))
{
yLower = y - 1
if (yLower == 0) yLower = max(world$Y)
yUpper = y + 1
if (yUpper > max(world$Y)) yUpper = 1
neighborValues = c()
neighborValues = c(neighborValues, world[world$X == xUpper & world$Y == y,]$Value)
neighborValues = c(neighborValues, world[world$X == xLower & world$Y == y,]$Value)
neighborValues = c(neighborValues, world[world$X == x & world$Y == yUpper,]$Value)
neighborValues = c(neighborValues, world[world$X == x & world$Y == yLower,]$Value)
if (includeDiagonals)
{
neighborValues = c(neighborValues, world[world$X == xUpper & world$Y == yUpper,]$Value)
neighborValues = c(neighborValues, world[world$X == xUpper & world$Y == yLower,]$Value)
neighborValues = c(neighborValues, world[world$X == xLower & world$Y == yUpper,]$Value)
neighborValues = c(neighborValues, world[world$X == xLower & world$Y == yLower,]$Value)
}
currentValue = world[newWorld$X == x & newWorld$Y == y,]$Value
newWorld[newWorld$X == x & newWorld$Y == y,]$Value = currentValue + updateFunction(currentValue, neighborValues)
}
}
newWorld$Density = newWorld$Value / sum(newWorld$Value)
return (newWorld)
}
# Update the world and display its new landscape
world = updateWorld(world, getAgentChange)
ggplot(world, aes(X, Y, z = Density)) +
geom_raster(aes(fill = Density)) +
geom_contour(colour = "white")
# Assembles the base data - what elevation is each coordinate at for each time step
buildTimeWorld <- function(width, height, iterations)
{
lastStep = buildWorld(width, height, seq(100, 100+(height * width * 0.1), 0.1))
lastStep$Time = rep(1, nrow(lastStep))
totalData = lastStep
for (t in 2:iterations)
{
newData = updateWorld(lastStep, getAgentChange)
newData$Time = rep(t, nrow(newData))
lastStep = newData
totalData = rbind(totalData, newData)
}
return(totalData)
}
# Given a set of points over time, outputs the optimal point at each time step
getOptimalData <- function(timeWorld)
{
optimalData = data.frame(OptimalX = 0, OptimalY = 0, Time = 0)
i = 1
for (t in min(timeWorld$Time):max(timeWorld$Time))
{
timeSubset = timeWorld[timeWorld$Time == t,]
optimalData[i,]$OptimalX = timeSubset[timeSubset$Value == max(timeSubset$Value),][1,]$X
optimalData[i,]$OptimalY = timeSubset[timeSubset$Value == max(timeSubset$Value),][1,]$Y
optimalData[i,]$Time = t
i = i + 1
}
return (optimalData)
}
# Given a set of points over time and an optimization function, returns the optimization function's
# guess of optimality for each time step
getGuesses <- function(timeWorld, optimizeFunction, ticksPerTime = 1, ...)
{
guessData = data.frame(GuessedX = 0, GuessedY = 0, Time = 0)
xGuess = floor(mean(timeWorld$X))
yGuess = floor(mean(timeWorld$Y))
i = 1
for (t in min(timeWorld$Time):max(timeWorld$Time))
{
timeSubset = timeWorld[timeWorld$Time == t,]
guesses = optimizeFunction(timeSubset, ticksPerTime, xGuess, yGuess, ...)
xGuess = guesses[1]
yGuess = guesses[2]
guessData[i,]$GuessedX = xGuess
guessData[i,]$GuessedY = yGuess
guessData[i,]$Time = t
i = i + 1
}
return (guessData)
}
# Horrible ugly monster function to encapsulate a bunch of work.
# Builds a world (or reuses the one passed in a timeWorld). Uses that to obtain optimal points
# (if includeOptimalPoints is true), as well as run an optimization function (if optimizeFunction is not
# null). Returns an object with three dataframes (or NULLs) - worldData, optimalPoints, and guessedPoints.
assembleData <- function (width = 17, height = 17, iterations = 50, timeWorld = NULL, optimizeFunction = NULL, ticksPerTime = 1, includeOptimalPoints = T, ...)
{
finalData <- list(worldData = NULL, optimalPoints = NULL, guessedPoints = NULL)
if (is.null(timeWorld))
{
finalData$worldData = buildTimeWorld(width, height, iterations)
} else {
finalData$worldData = timeWorld
}
if (!is.null(optimizeFunction))
{
finalData$guessedPoints = getGuesses(finalData$worldData, optimizeFunction, ticksPerTime, ...)
}
if (includeOptimalPoints)
{
finalData$optimalPoints = getOptimalData(finalData$worldData)
}
return (finalData)
}
# Assembling a world only
finalData = assembleData(16, 16, 50)
# Save out the world for reuse
timeWorld = finalData$worldData
# See how much faster this is now?
finalData = assembleData(timeWorld = timeWorld)
visualizeFunction <- function(allData, fromtime = 25, untiltime = 50)
{
fromtime = max(min(fromtime, max(allData$worldData$Time)), min(allData$worldData$Time))
untiltime = max(min(untiltime, max(allData$worldData$Time)), min(allData$worldData$Time))
worldData = allData$worldData[allData$worldData$Time >= fromtime & allData$worldData$Time <= untiltime,]
# Assemble base animation
animation = ggplot(worldData, aes(X, Y, z = Density)) +
geom_raster(aes(fill = Density)) +
transition_states(Time, transition_length = 10, state_length = 0, wrap=F)
labs(title = "Time: {closest_state}") +
enter_fade() +
exit_fade() +
ease_aes("linear")
# Add optimal points, if present
if (!is.null(allData$optimalPoints))
{
optData = allData$optimalPoints[allData$optimalPoints$Time >= fromtime & allData$optimalPoints$Time <= untiltime,]
lengthenFactor = nrow(worldData) / nrow(optData)
final = optData
for(i in 2:lengthenFactor) final = rbind(final, optData)
optData = final
optData = optData[with(optData, order(Time)), ]
animation = animation + geom_point(x=optData$OptimalX, y=optData$OptimalY, colour="hotpink", size=4)
}
# Add guessed points, if present
if (!is.null(allData$guessedPoints))
{
guessData = allData$guessedPoints[allData$guessedPoints$Time >= fromtime & allData$guessedPoints$Time <= untiltime,]
lengthenFactor = nrow(worldData) / nrow(guessData)
final = guessData
for(i in 2:lengthenFactor) final = rbind(final, guessData)
guessData = final
guessData = guessData[with(guessData, order(Time)), ]
animation = animation + geom_point(x=guessData$GuessedX, y=guessData$GuessedY, colour="lawngreen", size=4)
}
# Display animation
animation
}
visualizeFunction(finalData)
All functions have a standardized signature:
evaluate <- function (allData, fromtime = 25, untiltime = 50)
{
if (is.null(allData$optimalPoints) || is.null(allData$guessedPoints))
{
stop("Incomplete data for evaluate; allData is missing either optimalPoints or guessedPoints")
}
expected = allData$optimalPoints
actual = allData$guessedPoints
fromtime = max(min(fromtime, max(expected$Time)), min(expected$Time))
untiltime = max(min(untiltime, max(expected$Time)), min(expected$Time))
xLength = max(allData$worldData$X) - min(allData$worldData$X)
yLength = max(allData$worldData$Y) - min(allData$worldData$Y)
averageDistance = 0
for (t in fromtime:untiltime)
{
xLower = min(expected[expected$Time == t,]$OptimalX, actual[actual$Time == t,]$GuessedX)
xUpper = max(expected[expected$Time == t,]$OptimalX, actual[actual$Time == t,]$GuessedX)
yLower = min(expected[expected$Time == t,]$OptimalY, actual[actual$Time == t,]$GuessedY)
yUpper = max(expected[expected$Time == t,]$OptimalY, actual[actual$Time == t,]$GuessedY)
baseDistance = sqrt((xUpper - xLower) ^ 2 + (yUpper - yLower) ^ 2)
xWrap = sqrt((xUpper - (xLower + xLength)) ^ 2 + (yUpper - yLower) ^ 2)
yWrap = sqrt((xUpper - xLower) ^ 2 + (yUpper - (yLower + yLength)) ^ 2)
bothWrap = sqrt((xUpper - (xLower + xLength)) ^ 2 + (yUpper - (yLower + yLength)) ^ 2)
averageDistance = averageDistance + min(baseDistance, xWrap, yWrap, bothWrap)
}
averageDistance = averageDistance / (untiltime - fromtime + 1)
return (averageDistance)
}
randomGuesses <- function(currentData, ticks, startX = 1, startY = 1)
{
bestX = round(runif(1, min = min(currentData$X), max = max(currentData$X)))
bestY = round(runif(1, min = min(currentData$Y), max = max(currentData$Y)))
return (c(bestX, bestY))
}
hillClimb <- function(currentData, ticks, startX = 1, startY = 1)
{
bestX = startX
bestY = startY
bestValue = currentData[currentData$X == bestX & currentData$Y == bestY,]$Value
for (i in 1:ticks)
{
xLower = bestX - 1
if (xLower == 0) xLower = max(currentData$X)
xUpper = bestX + 1
if (xUpper > max(currentData$X)) xUpper = 1
yLower = bestY - 1
if (yLower == 0) yLower = max(currentData$Y)
yUpper = bestY + 1
if (yUpper > max(currentData$Y)) yUpper = 1
if (currentData[currentData$X == xLower & currentData$Y == bestY,]$Value > bestValue)
{
bestX = xLower
bestValue = currentData[currentData$X == xLower & currentData$Y == bestY,]$Value
}
if (currentData[currentData$X == xUpper & currentData$Y == bestY,]$Value > bestValue)
{
bestX = xUpper
bestValue = currentData[currentData$X == xUpper & currentData$Y == bestY,]$Value
}
if (currentData[currentData$X == bestX & currentData$Y == yUpper,]$Value > bestValue)
{
bestY = yUpper
bestValue = currentData[currentData$X == bestX & currentData$Y == yUpper,]$Value
}
if (currentData[currentData$X == bestX & currentData$Y == yLower,]$Value > bestValue)
{
bestY = yLower
bestValue = currentData[currentData$X == bestX & currentData$Y == yLower,]$Value
}
}
return (c(bestX, bestY))
}
Accepts two additional parameters:
linearCoolingFunction <- function (x) { return(1 - 0.1 * x) }
powerCoolingFunction <- function (x) { return(0.75 ^ x) }
decayCoolingFunction <- function (x) { return(1 / (x + 1)) }
simulatedAnnealing <- function(currentData, ticks, startX = 1, startY = 1, coolingFunction = NULL, temperatureTickAmount = NULL)
{
if (is.null(coolingFunction))
{
stop("Simulated Annealing requires CoolingFunction!")
}
if (is.null(temperatureTickAmount))
{
stop("Simulated Annealing requires TemperatureTickAmount!")
}
bestX = startX
bestY = startY
bestValue = currentData[currentData$X == bestX & currentData$Y == bestY,]$Value
currentCoolingValue = 0
for (i in 1:ticks)
{
xLower = bestX - 1
if (xLower == 0) xLower = max(currentData$X)
xUpper = bestX + 1
if (xUpper > max(currentData$X)) xUpper = 1
yLower = bestY - 1
if (yLower == 0) yLower = max(currentData$Y)
yUpper = bestY + 1
if (yUpper > max(currentData$Y)) yUpper = 1
nextX = bestX
nextY = bestY
nextValue = 0
if (currentData[currentData$X == xLower & currentData$Y == bestY,]$Value > nextValue)
{
nextX = xLower
nextValue = currentData[currentData$X == xLower & currentData$Y == bestY,]$Value
}
if (currentData[currentData$X == xUpper & currentData$Y == bestY,]$Value > nextValue)
{
nextX = xUpper
nextValue = currentData[currentData$X == xUpper & currentData$Y == bestY,]$Value
}
if (currentData[currentData$X == bestX & currentData$Y == yUpper,]$Value > nextValue)
{
nextY = yUpper
nextValue = currentData[currentData$X == bestX & currentData$Y == yUpper,]$Value
}
if (currentData[currentData$X == bestX & currentData$Y == yLower,]$Value > nextValue)
{
nextY = yLower
nextValue = currentData[currentData$X == bestX & currentData$Y == yLower,]$Value
}
# The part that makes this simulated annealing, not hill climbing
currentCoolingValue = currentCoolingValue + temperatureTickAmount
if ((nextValue < bestValue && runif(1) < coolingFunction(currentCoolingValue)) ||
nextValue > bestValue)
{
bestValue = nextValue
bestX = nextX
bestY = nextY
}
}
return (c(bestX, bestY))
}
randomOutput = assembleData(timeWorld = timeWorld, optimizeFunction = randomGuesses, ticksPerTime = 1)
cat("Random guesses were, on average, off by", evaluate(randomOutput), "units.\n")
## Random guesses were, on average, off by 5.882587 units.
visualizeFunction(randomOutput)
hillClimbOutput = assembleData(timeWorld = timeWorld, optimizeFunction = hillClimb, ticksPerTime = 10)
cat("Hill climbing was, on average, off by", evaluate(hillClimbOutput), "units.\n")
## Hill climbing was, on average, off by 6.64922 units.
visualizeFunction(hillClimbOutput)